6. Non-Cartesian spiral under-sampling

In this notebook, you can play with the design parameters to regenerate different spiral in-out patterns (so, we draw as many spiral arches as the number of shots). You can play with the number of shots by changing the under-sampling factor.

• Author: Philippe Ciuciu (philippe.ciuciu@cea.fr)
• Date: 06/24/2022
• Update: 01/29/2024 (dependence on MRI-NUFFT for density compensation)
• Update: 02/04/2025 (remove dependence on pysap + add illustration on T2w contrast)
• Target: IEEE EMBS-SPS Summer School on Novel acquisition and image reconstruction strategies in accelerated Magnetic Resonance Imaging

#DISPLAY T2* MR IMAGE
'''
#DISPLAY BRAIN PHANTOM
在医学成像（如 MRI、CT）和信号处理中，幻影（Phantom）是指用于测试、校准或研究的人工数据或物理模型。
它通常是一个合成（模拟）图像或实验装置，用来模拟真实世界的结构，比如大脑、心脏或其他人体组织。
'''
# 1. 导入必要的库
%matplotlib inline

import numpy as np
import os.path as op
import os
import math ; import cmath
import matplotlib.pyplot as plt
import sys

import brainweb_dl as bwdl                                          # BrainWeb MRI 数据下载
import mrinufft                                                     # MRI 非均匀 FFT (NUFFT) 处理库
from mrinufft.trajectories import display
from mrinufft import get_density                                    # 计算密度补偿
from mrinufft.density import voronoi
from mrinufft.trajectories import initialize_2D_spiral
from mrinufft.trajectories.display import display_2D_trajectory     # MRI 轨迹可视化


from skimage import data, img_as_float, io, filters                 # 处理 MRI 图像
from modopt.math.metrics import ssim                                # 计算结构相似度 SSIM

# 2. 设置 Matplotlib 显示参数
plt.rcParams["image.origin"]="lower"                                # 设置图像的原点在左下角（默认是左上角）。
plt.rcParams["image.cmap"]='Greys_r'                                # 设置图像颜色映射为灰度反转（黑白）。

# 3. 下载并处理 MRI 图像
mri_img = bwdl.get_mri(4, "T1")[70, ...].astype(np.float32)
'''
bwdl.get_mri(4, "T1") : 从 BrainWeb 数据库下载一个 T1 加权 MRI 扫描
70 : 获取 MRI 体数据（3D）中的第 70 层切片，即 2D 切片图像。
.astype(np.float32) : 转换数据类型为 float32，以便进行计算。
''' 

# 4. 显示 MRI 图像
print(mri_img.shape)                                                # (256, 256)
img_size = mri_img.shape[0]                                         # 获取图像的宽度（假设是正方形图像）。

# 5. 绘制 MRI 图像
plt.figure()
plt.imshow(abs(mri_img))                                            # 显示 MRI 图像，abs() 主要用于防止负像素值影响显示。
plt.title("Original brain image")
plt.show()

(256, 256)

Code given by Github

# set up the first shot
rfactor = 4
num_shots = math.ceil(img_size/rfactor)
print("number of shots: {}".format(num_shots))

# define the regularly spaced samples on a single shot
#nsamples = (np.arange(0,img_size) - img_size//2)/(img_size)
num_samples = img_size
num_samples = (num_samples + 1) // 2
print("number of samples: {}".format(num_samples))
num_revolutions = 1

shot = np.arange(0, num_samples, dtype=np.complex_)
radius = shot / num_samples * 1 / (2 * np.pi) * (1 - np.finfo(float).eps)
angle = np.exp(2 * 1j * np.pi * shot / num_samples * num_revolutions)
# first half of the spiral
single_shot = np.multiply(radius, angle)
# add second half of the spiral
#single_shot = np.append(np.flip(single_shot, axis=0), -single_shot[1:])
single_shot = np.append(np.flip(single_shot, axis=0), -single_shot)
print("number of samples per shot: {}".format(np.size(single_shot)))

k_shots = np.array([], dtype = np.complex_)
#for i in vec_shots:
for i in np.arange(0, num_shots):
    shot_rotated = single_shot * np.exp(1j * 2 * np.pi * i / (num_shots * 2))
    k_shots = np.append(k_shots, shot_rotated)

print(k_shots.shape)
kspace_loc = np.zeros((len(k_shots),2))
kspace_loc[:,0] = k_shots.real
kspace_loc[:,1] = k_shots.imag

#Plot full initialization
kspace = plt.figure(figsize = (8,8))
#plot shots
plt.scatter(kspace_loc[::4,0], kspace_loc[::4,1], marker = '.')
plt.title("Spiral undersampling R = %d" %rfactor)

axes = plt.gca()
plt.grid()

number of shots: 64
number of samples: 128
number of samples per shot: 256
(16384,)

print(np.arange(0, num_shots))

[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63]

import warnings
warnings.filterwarnings("ignore")  # 忽略所有警告

#data=convert_locations_to_mask(kspace_loc, mri_img.shape)

NufftOperator = mrinufft.get_operator("finufft")

# For better image quality we use a density compensation
density = "voronoi"
#density = "cell_count"
#ensity = None
# And create the associated operator.
if density == None:
    nufft = NufftOperator(
        kspace_loc.reshape(-1, 2), shape=mri_img.shape, density=None, n_coils=1
    )
elif density == "voronoi":
    voronoi_weights = get_density("voronoi", kspace_loc.reshape(-1, 2))
    nufft = NufftOperator(
        kspace_loc.reshape(-1, 2), shape=mri_img.shape, density=voronoi_weights, n_coils=1
    )
else:
    cell_count_weights = get_density("cell_count", kspace_loc.reshape(-1, 2), shape=mri_img.shape, osf=2.0)
    nufft = NufftOperator(
        kspace_loc.reshape(-1, 2), shape=mri_img.shape, density=cell_count_weights, n_coils=1, upsampfac=2.0
    )

kspace_obs = nufft.op(mri_img)  # Image -> Kspace
image2 = nufft.adj_op(kspace_obs)  # Kspace -> Image

grid_space = np.linspace(-0.5, 0.5, num=mri_img.shape[0])
grid2D = np.meshgrid(grid_space, grid_space)
#rid_soln = gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs,
#                                                tuple(grid2D), 'linear')
plt.imshow(np.abs(image2), cmap='gray')
# Calculate SSIM
base_ssim = ssim(image2, mri_img)
plt.title('Gridded Solution\nSSIM = ' + str(base_ssim))
plt.show()

Q1  

• Second block: Vary parameter rfactor (from 4 to 32) in the definition of the number of shots (i.e. spokes) to perform spiral under-sampling and observe its impact on image quality (SSIM score).

def Reconstruction_v1 (mask, rfactor, density):
    '''
    mask = "Spiral"
    rfactor : from 4 to 32
    #density = "voronoi"
    #density = "cell_count"
    #density = "None"
    '''
    if mask == "Radial" :
        # set up the first shot
        ## rfactor = 4
        nb_shots = math.ceil(img_size/rfactor)
        ## print("number of shots: {}".format(nb_shots))

        # vectorize the nb of shots
        vec_shots = np.arange(0,nb_shots)

        # define the regularly spaced samples on a single shot
        nsamples = (np.arange(0,img_size) - img_size//2)/(img_size)
        ## print("number of samples per shot: {}".format(np.size(nsamples)))

        # NumPy 2.0 : np.complex_ -> np.complex128
        shot_c = np.array(nsamples, dtype=np.complex128)
        shots = np.array([], dtype=np.complex128)
        #shot_c = np.array(nsamples, dtype = np.complex_)
        #shots = np.array([], dtype = np.complex_)
        # acculumate shots after rotating the initial one by the right angular increment
        for k in vec_shots:
            shots = np.append(shots, shot_c * np.exp(2 * np.pi * 1j * k/(2*nb_shots)))

        kspace_loc = np.zeros((len(shots),2))

        #assign real and imaginary parts of complex-valued k-space trajectories to k-space locations
        kspace_loc[:,0] = shots.real
        kspace_loc[:,1] = shots.imag

        '''
        #Plot full initialization
        kspace = plt.figure(figsize = (8,8))
        ##plot shots
        plt.scatter(kspace_loc[:,0],kspace_loc[:,1], marker = '.')
        plt.title("Radial undersampling R = %d" %rfactor)

        axes = plt.gca()
        plt.grid()
        '''

        #kspace_loc = mrinufft.initialize_2D_radial(Nc=nb_shots, Ns=nsamples)

        ## The real deal starts here ##
        # Choose your NUFFT backend (installed independly from the package)
        NufftOperator = mrinufft.get_operator("finufft")

        # For better image quality we use a density compensation
        #density = "voronoi"
        #density = "cell_count"
        #density = None
        # And create the associated operator.
        if density == None:
            nufft = NufftOperator(
                kspace_loc.reshape(-1, 2), shape=mri_img.shape, density=None, n_coils=1
            )
        elif density == "voronoi":
            voronoi_weights = get_density("voronoi", kspace_loc.reshape(-1, 2))
            nufft = NufftOperator(
                kspace_loc.reshape(-1, 2), shape=mri_img.shape, density=voronoi_weights, n_coils=1
            )
        else:
            cell_count_weights = get_density("cell_count", kspace_loc.reshape(-1, 2), shape=mri_img.shape, osf=2.0)
            nufft = NufftOperator(
                kspace_loc.reshape(-1, 2), shape=mri_img.shape, density=cell_count_weights, n_coils=1, upsampfac=2.0
            )

        kspace_obs = nufft.op(mri_img)  # Image -> Kspace
        image2 = nufft.adj_op(kspace_obs)  # Kspace -> Image

        '''
        grid_space = np.linspace(-0.5, 0.5, num=mri_img.shape[0])
        grid2D = np.meshgrid(grid_space, grid_space)
        #grid_soln = gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs,
        #                                                 tuple(grid2D), 'linear')
        plt.imshow(np.abs(image2), cmap='gray')
        # Calculate SSIM
        base_ssim = ssim(image2, mri_img)
        plt.title('Adjoint NUFFT solution\nSSIM = ' + str(base_ssim))
        plt.show()
        '''
        # Calculate SSIM
        base_ssim = ssim(image2, mri_img)
        
    elif mask == "Spiral" :
        # set up the first shot
        ## rfactor = 4
        num_shots = math.ceil(img_size/rfactor)
        ## print("number of shots: {}".format(num_shots))

        # define the regularly spaced samples on a single shot
        #nsamples = (np.arange(0,img_size) - img_size//2)/(img_size)
        num_samples = img_size
        num_samples = (num_samples + 1) // 2
        ## print("number of samples: {}".format(num_samples))
        num_revolutions = 1

        shot = np.arange(0, num_samples, dtype=np.complex_)
        radius = shot / num_samples * 1 / (2 * np.pi) * (1 - np.finfo(float).eps)
        angle = np.exp(2 * 1j * np.pi * shot / num_samples * num_revolutions)
        # first half of the spiral
        single_shot = np.multiply(radius, angle)
        # add second half of the spiral
        #single_shot = np.append(np.flip(single_shot, axis=0), -single_shot[1:])
        single_shot = np.append(np.flip(single_shot, axis=0), -single_shot)
        ## print("number of samples per shot: {}".format(np.size(single_shot)))

        k_shots = np.array([], dtype = np.complex_)
        #for i in vec_shots:
        for i in np.arange(0, num_shots):
            shot_rotated = single_shot * np.exp(1j * 2 * np.pi * i / (num_shots * 2))
            k_shots = np.append(k_shots, shot_rotated)

        ## print(k_shots.shape)
        kspace_loc = np.zeros((len(k_shots),2))
        kspace_loc[:,0] = k_shots.real
        kspace_loc[:,1] = k_shots.imag

        '''
        #Plot full initialization
        kspace = plt.figure(figsize = (8,8))
        #plot shots
        plt.scatter(kspace_loc[::4,0], kspace_loc[::4,1], marker = '.')
        plt.title("Spiral undersampling R = %d" %rfactor)

        axes = plt.gca()
        plt.grid()
        '''
        
        #data=convert_locations_to_mask(kspace_loc, mri_img.shape)

        NufftOperator = mrinufft.get_operator("finufft")

        # For better image quality we use a density compensation
        density = "voronoi"
        #density = "cell_count"
        #ensity = None
        # And create the associated operator.
        if density == None:
            nufft = NufftOperator(
                kspace_loc.reshape(-1, 2), shape=mri_img.shape, density=None, n_coils=1
            )
        elif density == "voronoi":
            voronoi_weights = get_density("voronoi", kspace_loc.reshape(-1, 2))
            nufft = NufftOperator(
                kspace_loc.reshape(-1, 2), shape=mri_img.shape, density=voronoi_weights, n_coils=1
            )
        else:
            cell_count_weights = get_density("cell_count", kspace_loc.reshape(-1, 2), shape=mri_img.shape, osf=2.0)
            nufft = NufftOperator(
                kspace_loc.reshape(-1, 2), shape=mri_img.shape, density=cell_count_weights, n_coils=1, upsampfac=2.0
            )

        kspace_obs = nufft.op(mri_img)  # Image -> Kspace
        image2 = nufft.adj_op(kspace_obs)  # Kspace -> Image

        '''
        grid_space = np.linspace(-0.5, 0.5, num=mri_img.shape[0])
        grid2D = np.meshgrid(grid_space, grid_space)
        #rid_soln = gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs,
        #                                                tuple(grid2D), 'linear')
        plt.imshow(np.abs(image2), cmap='gray')
        # Calculate SSIM
        base_ssim = ssim(image2, mri_img)
        plt.title('Gridded Solution\nSSIM = ' + str(base_ssim))
        plt.show()
        '''
        # Calculate SSIM
        base_ssim = ssim(image2, mri_img)
        
    return kspace_loc, image2, base_ssim

$\leadsto$   Here we write a function using a mask based on Spiral Undersampling in k-space and reconstructs the MRI image using the Non-Uniform Fast Fourier Transform (NUFFT). The undersampling is controlled by the undersampling factor rfactor and density compensation method density.

It generates Radial k-space Sampling Locations : It defines total number of radial shots in k-space nb_shots, also the sample points along a single radial line nsamples normalized between [-0.5, 0.5]; We generate a single spiral arm : The radius is computed as radius = shot / num_samples * 1 / (2 * np.pi) * (1 - np.finfo(float).eps), The angle is defined as angle = np.exp(2 * 1j * np.pi * shot / num_samples * num_revolutions), We generate the first half of the spiral and mirror it to ensure symmetry; We rotate each radial line by an increment of 2 * np.pi / (2 * nb_shots) ; We store the real kspace_loc[:,0] and imaginary kspace_loc[:,1] components of k-space locations.

$\leadsto$   NUFFT adapts to non-uniformly spaced data by modifying the basis functions :

$$ F_{\text{NUFFT}} = DWF $$

where:

• $D$ is the density compensation

• $W$ is the interpolation function

• $F$ is the standard Fourier transform, or

  $$ F(u,v) = \sum_{x,y} f(x,y)e^{-2\pi i (ux+vy)} $$

We choose the NUFFT (Non-Uniform FFT) backend $W$ : NufftOperator = mrinufft.get_operator("finufft") ; We apply density compensation $D$ ("voronoi" or "cell_count") if needed.

We compute non-uniform Fourier transform of the MRI image $W$ : kspace_obs = nufft.op(mri_img) ; And we reconstruct the image from the undersampled k-space data $W$ by image2 = nufft.adj_op(kspace_obs)

mask="Spiral"
vector_rfactor = np.linspace(4, 32, 8)
vector_density = ["voronoi", "cell_count", "None"]

fig, axs = plt.subplots(len(vector_rfactor), 2, figsize=(24, 24))
for i in range(len(vector_rfactor)) :
    for j in range(2):

        mask = "Spiral"
        rfactor = vector_rfactor[i]
        density = vector_density[0]

        kspace_loc, image2, base_ssim = Reconstruction_v1(mask, rfactor, density)

        # Plot full initialization
        axs[i, 0].scatter(kspace_loc[:,0],kspace_loc[:,1], marker = '.')
        axs[i, 0].set_aspect('equal')  # 设置轴比例相等，保持正方形

        # 显示重建的图像
        axs[i, 1].imshow(np.abs(image2), cmap='Greys_r')

        axs[i, 0].set_xlabel(f"R = {rfactor}")
        # axs[i, 1].set_xlabel(f"SSIM = {base_ssim}")
        axs[i, 1].set_xlabel(f"SSIM = {base_ssim:.2f}")  # 保留2位小数
        axs[i, j].xaxis.set_label_position('top')  # 移动 xlabel 到顶部
        axs[i, j].set_ylabel(f"density = {density}")

axs[0, 0].set_title("Spiral undersampling")
axs[0, 1].set_title("Adjoint NUFFT solution")

# Add a main title (suptitle) to the figure
fig.suptitle(f"k-space locations and its Reconstruction Results for Mask: {mask}", fontsize=16)

# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96])  # Adjust layout to fit suptitle
plt.show()

$\leadsto$   In the case where mask = Radial, density = "voronoi", we plot 8 reconstructed images with rfactor varying in vector_rfactor = np.linspace(4, 32, 8) :

• Initially, as rfactor increases, the quality of the reconstructed image improves until rfactor ≃ 6, after which the image quality starts to degrade.

• To quantitatively evaluate the reconstruction accuracy, we can calculate the Structural Similarity Index between the true image and the reconstructed image.

$\leadsto$   We can compute the Structural Similarity Index to evaluate the quality of reconstruction :

$$ SSIM(x,y) = \frac{(2m_x m_y + C_1)(2\Sigma_{xy} + C_2)}{(m_x^2 + m_y^2 + C_1)(\Sigma_{xx}^2 + \Sigma_{yy}^2 + C_2)} $$

where $m$ represents the mean, and $\Sigma$ denotes the covariance, and $C_1$ and $C_2$ are constants.

The result represents the overall similarity between the two images, ranging from [0, 1], where 1 indicates they are identical, and 0 means there is no similarity at all.

In Python, import the module like :

from modopt.math.metrics import ssim

Then, simply use the following line to compute the SSIM score and map :

base_ssim = ssim(image2, mri_img)

mask="Spiral"
vector_rfactor = np.linspace(4, 32, 29)
vector_density = ["voronoi", "cell_count", "None"]

vector_base_ssim = []
vector_image2 = []

for i in range(len(vector_rfactor)) :

        mask = "Spiral"
        rfactor = vector_rfactor[i]
        density = vector_density[0]

        kspace_loc, image2, base_ssim = Reconstruction_v1(mask, rfactor, density)

        vector_image2.append(image2)
        vector_base_ssim.append(base_ssim)

plt.figure(figsize=(10, 6))
plt.plot(vector_rfactor , vector_base_ssim, linestyle='-', color='red', marker='.', markersize=8, markerfacecolor='red', label='Base SSIM')

plt.xlabel('Rfactor')
plt.ylabel('Base SSIM')
plt.title('Base SSIM vs Rfactor')
plt.grid(True)
plt.legend()
plt.show()

The results show that the maximum Base SSIM occurs at rfactor = 6 . Below, we display the reconstructed images for it :

• rfactor = 6 (corresponding to the maximum Base SSIM)

with density = "voronoi".

# Find the index and value of the ；maximum Frobenius norm.
max_index = np.argmax(vector_base_ssim)
image_corresponding_max_index = vector_image2[max_index]
rfactor_corresponding_max_index = vector_rfactor[max_index]
# print(rfactor_corresponding_max_index)

fig, axes = plt.subplots(1, 2, figsize=(8, 4))

axes[0].imshow(np.abs(image_corresponding_max_index), cmap='Greys_r')
axes[0].set_title(f'Rfactor (for Max Base SSIM) : {rfactor_corresponding_max_index:.2f}')
axes[0].axis('off')

axes[1].imshow(mri_img, cmap='Greys_r')
axes[1].set_title("True image")
axes[1].axis('off')

plt.tight_layout()
plt.show()

from scipy.spatial import Voronoi, voronoi_plot_2d

vor = Voronoi(kspace_loc.reshape(-1, 2))
fig = voronoi_plot_2d(vor)

fig = voronoi_plot_2d(vor, show_vertices=False, line_colors='orange',
                      line_width=2, line_alpha=0.6, point_size=2)
plt.show()

mri_2D = bwdl.get_mri(4, "T2")[150, ...].astype(np.float32)
print(mri_2D.shape)
plt.figure()
plt.imshow(abs(mri_2D))
plt.show()

(434, 362)

from mrinufft import get_density, get_operator, check_backend

traj = initialize_2D_spiral(256, 256).astype(np.float32)

nufft = get_operator("finufft")(traj, mri_2D.shape, density=False)
kspace = nufft.op(mri_2D)
adjoint = nufft.adj_op(kspace)

fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(abs(mri_2D))
display_2D_trajectory(traj, subfigure=axs[1])
axs[2].imshow(abs(adjoint))

voronoi_weights = get_density("voronoi", traj)
nufft_voronoi = get_operator("finufft")(
    traj, shape=mri_2D.shape, density=voronoi_weights
)
adjoint_voronoi = nufft_voronoi.adj_op(kspace)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(abs(mri_2D))
display_2D_trajectory(traj, subfigure=axs[1])
axs[2].imshow(abs(adjoint_voronoi))

flat_traj = traj.reshape(-1, 2)
weights = np.sqrt(np.sum(flat_traj**2, axis=1))
nufft = get_operator("finufft")(traj, shape=mri_2D.shape, density=weights)
adjoint_manual = nufft.adj_op(kspace)
fig, axs = plt.subplots(1, 3, figsize=(15, 5))
axs[0].imshow(abs(mri_2D))
axs[0].set_title("Ground Truth")
axs[1].imshow(abs(adjoint))
axs[1].set_title("no density compensation")
axs[2].imshow(abs(adjoint_manual))
axs[2].set_title("manual density compensation")

Text(0.5, 1.0, 'manual density compensation')

Q2  

• Activate the density compensation mechanism using the Voronoi diagram or cell counting technique and analyze how image quality evolves accordingly.

$\leadsto$   No density compensation mechanism : nufft = get_operator("finufft")(traj, mri_2D.shape, density=False)

TThe reconstructed image exhibits saturation and blurring. The quality of the reconstruction is poor.

$\leadsto$   When we choose Voronoi Diagram Based Density Compensation : voronoi_weights = get_density("voronoi", traj)

The reconstruction quality is significantly improved compared to the case with no density compensation.

$\leadsto$   When we choose manual density compensation flat_traj = traj.reshape(-1, 2) weights = np.sqrt(np.sum(flat_traj**2, axis=1))

The reconstruction quality is enhanced similarly to the Voronoi Diagram Based Density Compensation approach.

Q3  

• For a given value of parameter rfactor comment the differences in terms of image quality between radial and spiral undersampling.

vector_mask = ["Radial", "Spiral"]
vector_rfactor = np.linspace(4, 32, 8)
vector_density = ["voronoi", "cell_count", "None"]

fig, axs = plt.subplots(len(vector_rfactor), 4, figsize=(24, 24))
for i in range(len(vector_rfactor)) :
    for j in [0,1] :

        mask = "Radial"
        rfactor = vector_rfactor[i]
        density = vector_density[0]

        kspace_loc, image2, base_ssim = Reconstruction_v1(mask, rfactor, density)

        # Plot full initialization
        axs[i, 0].scatter(kspace_loc[:,0],kspace_loc[:,1], marker = '.')
        axs[i, 0].set_aspect('equal')  # 设置轴比例相等，保持正方形

        # 显示重建的图像
        axs[i, 1].imshow(np.abs(image2), cmap='Greys_r')

        axs[i, 0].set_xlabel(f"R = {rfactor}")
        # axs[i, 1].set_xlabel(f"SSIM = {base_ssim}")
        axs[i, 1].set_xlabel(f"SSIM = {base_ssim:.2f}")  # 保留2位小数
        axs[i, j].xaxis.set_label_position('top')  # 移动 xlabel 到顶部
        axs[i, j].set_ylabel(f"density = {density}")

    for j in [2,3] :

        mask = "Spiral"
        rfactor = vector_rfactor[i]
        density = vector_density[0]

        kspace_loc, image2, base_ssim = Reconstruction_v1(mask, rfactor, density)

        # Plot full initialization
        axs[i, 2].scatter(kspace_loc[:,0],kspace_loc[:,1], marker = '.')
        axs[i, 2].set_aspect('equal')  # 设置轴比例相等，保持正方形

        # 显示重建的图像
        axs[i, 3].imshow(np.abs(image2), cmap='Greys_r')

        axs[i, 2].set_xlabel(f"R = {rfactor}")
        # axs[i, 1].set_xlabel(f"SSIM = {base_ssim}")
        axs[i, 3].set_xlabel(f"SSIM = {base_ssim:.2f}")  # 保留2位小数
        axs[i, j].xaxis.set_label_position('top')  # 移动 xlabel 到顶部
        axs[i, j].set_ylabel(f"density = {density}")

axs[0, 0].set_title("Radial undersampling")
axs[0, 1].set_title("Adjoint NUFFT solution")

axs[0, 2].set_title("Spiral undersampling")
axs[0, 3].set_title("Adjoint NUFFT solution")

# Add a main title (suptitle) to the figure
fig.suptitle(f"k-space locations and its Reconstruction Results for Mask: {vector_mask[0]} and {vector_mask[1]}", fontsize=16)

# Display the figure with adjusted layout
plt.tight_layout(rect=[0, 0, 1, 0.96])  # Adjust layout to fit suptitle
plt.show()

$\leadsto$   For the same undersampling factor R, Spiral underampling yields a higher SSIM (or the quality of reconstruction is better) compared to Radial undersampling .

This may be due to the shape of the undersampling in k-space: Spiral undersampling provides greater frequency variability, while Radial undersampling consists of straight lines, contains less information, and may include some symmetric components.